Working with SWOT Level 2 Water Mask Raster Image Data Product:

Working with SWOT Level 2 Water Mask Raster Image Data Product:

In AWS Cloud Version

Author: Nicholas Tarpinian, PO.DAAC

Summary & Learning Objectives

Notebook showcasing how to work with SWOT Raster Image data product version 2.0 on a local machine

  • Utilizing the earthaccess Python package. For more information visit: https://nsidc.github.io/earthaccess/
  • Option to query the new dataset based on users choice; choosing between two resolutions either by ‘100m’ or ‘250m’.
  • Visualizing the dataset by choosing one of its various variables.
  • Visualizing multiple raster images on a single map.
  • Stacking multiple raster images and creating a time dimenstion to analyze over time.

Requirements

1. Compute environment

This tutorial is written to run in the following environment: - AWS instance running in us-west-2: NASA Earthdata Cloud data in S3 can be directly accessed via and s3fs session; this access is limited to requests made within the US West (Oregon) (code: us-west-2) AWS region.

2. Earthdata Login

An Earthdata Login account is required to access data, as well as discover restricted data, from the NASA Earthdata system. Thus, to access NASA data, you need Earthdata Login. Please visit https://urs.earthdata.nasa.gov to register and manage your Earthdata Login account. This account is free to create and only takes a moment to set up.

Import libraries

%matplotlib inline

import os
import io
import s3fs
import xarray as xr
import numpy as np
from datetime import datetime
from os.path import isfile, basename, abspath
from pathlib import Path
import matplotlib.pyplot as plt
from matplotlib import animation
from matplotlib.animation import FuncAnimation, PillowWriter
from IPython.display import display, Image, HTML
import hvplot
import hvplot.xarray 
import holoviews as hv
from holoviews import opts
import earthaccess
from earthaccess import Auth, DataCollections, DataGranules, Store

hv.extension('bokeh')

Authentication with earthaccess

In this notebook, we will be calling the authentication in the below cell.

auth = earthaccess.login(strategy="interactive", persist=True)

Search using earthaccess for SWOT Raster products

Each dataset has it’s own unique collection concept ID. For the SWOT_L2_HR_Raster_2.0 dataset, we can find the collection ID here.

For this tutorial, we are looking at the Lake Mead Reservoir.

We used bbox finder to get the exact coordinates for our area of interest.

raster_results = earthaccess.search_data(
    short_name = 'SWOT_L2_HR_RASTER_2.0',
    bounding_box=(-115.112686,35.740939,-114.224167,36.937819),
    temporal =('2024-01-01 00:00:00', '2024-01-05 23:59:59'),
    granule_name = '*_100m_*',
    count =200
)
Granules found: 2
fs_s3 = earthaccess.get_s3fs_session(results=raster_results)

Visualizing the Dataset

Let’s now visualize an individual layer for a single file path using Xarray to read the NetCDF image.

raster_link = earthaccess.results.DataGranule.data_links(raster_results[1], access='direct')[0]
s3_file_obj5 = fs_s3.open(raster_link, mode='rb')
ds_raster = xr.open_dataset(s3_file_obj5, engine='h5netcdf')
ds_raster
<xarray.Dataset>
Dimensions:                  (x: 1505, y: 1506)
Coordinates:
  * x                        (x) float64 5.931e+05 5.932e+05 ... 7.435e+05
  * y                        (y) float64 4.026e+06 4.026e+06 ... 4.176e+06
Data variables: (12/39)
    crs                      object ...
    longitude                (y, x) float64 ...
    latitude                 (y, x) float64 ...
    wse                      (y, x) float32 ...
    wse_qual                 (y, x) float32 ...
    wse_qual_bitwise         (y, x) float64 ...
    ...                       ...
    load_tide_fes            (y, x) float32 ...
    load_tide_got            (y, x) float32 ...
    pole_tide                (y, x) float32 ...
    model_dry_tropo_cor      (y, x) float32 ...
    model_wet_tropo_cor      (y, x) float32 ...
    iono_cor_gim_ka          (y, x) float32 ...
Attributes: (12/49)
    Conventions:                   CF-1.7
    title:                         Level 2 KaRIn High Rate Raster Data Product
    source:                        Ka-band radar interferometer
    history:                       2024-01-05T10:46:43Z : Creation
    platform:                      SWOT
    references:                    V1.1.1
    ...                            ...
    x_min:                         593100.0
    x_max:                         743500.0
    y_min:                         4025600.0
    y_max:                         4176100.0
    institution:                   CNES
    product_version:               01
ds_raster.wse.hvplot.image(y='y', x='x').opts(width=1000, height=600, colorbar=True)

Visualizing Multiple Images

Let’s now visualize the multiple raster images that were searched and explore the data.

Utilizing xarray.open_mfdataset which supports the opening of multiple files.

ds = xr.open_mfdataset(earthaccess.open(raster_results), engine='h5netcdf',combine='nested', concat_dim='x')
ds
Opening 2 granules, approx size: 0.06 GB
using endpoint: https://archive.swot.podaac.earthdata.nasa.gov/s3credentials
<xarray.Dataset>
Dimensions:                  (y: 2760, x: 3009)
Coordinates:
  * y                        (y) float64 3.9e+06 3.9e+06 ... 4.176e+06 4.176e+06
  * x                        (x) float64 5.683e+05 5.684e+05 ... 7.435e+05
Data variables: (12/39)
    crs                      (x) object b'1' b'1' b'1' b'1' ... b'1' b'1' b'1'
    longitude                (y, x) float64 dask.array<chunksize=(502, 502), meta=np.ndarray>
    latitude                 (y, x) float64 dask.array<chunksize=(502, 502), meta=np.ndarray>
    wse                      (y, x) float32 dask.array<chunksize=(752, 752), meta=np.ndarray>
    wse_qual                 (y, x) float32 dask.array<chunksize=(2760, 1504), meta=np.ndarray>
    wse_qual_bitwise         (y, x) float64 dask.array<chunksize=(752, 752), meta=np.ndarray>
    ...                       ...
    load_tide_fes            (y, x) float32 dask.array<chunksize=(752, 752), meta=np.ndarray>
    load_tide_got            (y, x) float32 dask.array<chunksize=(752, 752), meta=np.ndarray>
    pole_tide                (y, x) float32 dask.array<chunksize=(752, 752), meta=np.ndarray>
    model_dry_tropo_cor      (y, x) float32 dask.array<chunksize=(752, 752), meta=np.ndarray>
    model_wet_tropo_cor      (y, x) float32 dask.array<chunksize=(752, 752), meta=np.ndarray>
    iono_cor_gim_ka          (y, x) float32 dask.array<chunksize=(752, 752), meta=np.ndarray>
Attributes: (12/49)
    Conventions:                   CF-1.7
    title:                         Level 2 KaRIn High Rate Raster Data Product
    source:                        Ka-band radar interferometer
    history:                       2024-01-05T11:01:19Z : Creation
    platform:                      SWOT
    references:                    V1.1.1
    ...                            ...
    x_min:                         568300.0
    x_max:                         718600.0
    y_min:                         3900200.0
    y_max:                         4050500.0
    institution:                   CNES
    product_version:               01
    • y
      PandasIndex
      PandasIndex(Float64Index([3900200.0, 3900300.0, 3900400.0, 3900500.0, 3900600.0, 3900700.0,
                    3900800.0, 3900900.0, 3901000.0, 3901100.0,
                    ...
                    4175200.0, 4175300.0, 4175400.0, 4175500.0, 4175600.0, 4175700.0,
                    4175800.0, 4175900.0, 4176000.0, 4176100.0],
                   dtype='float64', name='y', length=2760))
    • x
      PandasIndex
      PandasIndex(Float64Index([568300.0, 568400.0, 568500.0, 568600.0, 568700.0, 568800.0,
                    568900.0, 569000.0, 569100.0, 569200.0,
                    ...
                    742600.0, 742700.0, 742800.0, 742900.0, 743000.0, 743100.0,
                    743200.0, 743300.0, 743400.0, 743500.0],
                   dtype='float64', name='x', length=3009))
  • Conventions :
    CF-1.7
    title :
    Level 2 KaRIn High Rate Raster Data Product
    source :
    Ka-band radar interferometer
    history :
    2024-01-05T11:01:19Z : Creation
    platform :
    SWOT
    references :
    V1.1.1
    reference_document :
    JPL D-56416 - Revision B - August 17, 2023
    contact :
    alexander.t.corben[at]jpl.nasa.gov
    cycle_number :
    8
    pass_number :
    511
    scene_number :
    109
    tile_numbers :
    [217 218 219 217 218 219]
    tile_names :
    511_217L, 511_218L, 511_219L, 511_217R, 511_218R, 511_219R
    tile_polarizations :
    V, V, V, H, H, H
    coordinate_reference_system :
    Universal Transverse Mercator
    resolution :
    100.0
    short_name :
    L2_HR_Raster
    descriptor_string :
    100m_UTM11S_N_x_x_x
    crid :
    PIC0
    pge_name :
    PGE_L2_HR_RASTER
    pge_version :
    5.0.4
    time_granule_start :
    2024-01-01T12:43:28.040906Z
    time_granule_end :
    2024-01-01T12:43:47.497154Z
    time_coverage_start :
    2024-01-01T12:43:28.041873Z
    time_coverage_end :
    2024-01-01T12:43:46.958982Z
    geospatial_lon_min :
    -116.24753568370588
    geospatial_lon_max :
    -114.56382182347716
    geospatial_lat_min :
    35.225727504637376
    geospatial_lat_max :
    36.59518831179102
    left_first_longitude :
    -116.24753568370588
    left_first_latitude :
    35.46460890403489
    left_last_longitude :
    -115.95968454013911
    left_last_latitude :
    36.59518831179102
    right_first_longitude :
    -114.86994488403649
    right_first_latitude :
    35.225727504637376
    right_last_longitude :
    -114.56382182347716
    right_last_latitude :
    36.35099280753312
    xref_l2_hr_pixc_files :
    SWOT_L2_HR_PIXC_008_511_217L_20240101T124328_20240101T124337_PIC0_01.nc, SWOT_L2_HR_PIXC_008_511_218L_20240101T124336_20240101T124347_PIC0_01.nc, SWOT_L2_HR_PIXC_008_511_219L_20240101T124346_20240101T124357_PIC0_01.nc, SWOT_L2_HR_PIXC_008_511_217R_20240101T124328_20240101T124337_PIC0_01.nc, SWOT_L2_HR_PIXC_008_511_218R_20240101T124336_20240101T124347_PIC0_01.nc, SWOT_L2_HR_PIXC_008_511_219R_20240101T124346_20240101T124357_PIC0_01.nc
    xref_l2_hr_pixcvec_files :
    SWOT_L2_HR_PIXCVec_008_511_217L_20240101T124328_20240101T124337_PIC0_01.nc, SWOT_L2_HR_PIXCVec_008_511_218L_20240101T124336_20240101T124347_PIC0_01.nc, SWOT_L2_HR_PIXCVec_008_511_219L_20240101T124346_20240101T124357_PIC0_01.nc, SWOT_L2_HR_PIXCVec_008_511_217R_20240101T124328_20240101T124337_PIC0_01.nc, SWOT_L2_HR_PIXCVec_008_511_218R_20240101T124336_20240101T124347_PIC0_01.nc, SWOT_L2_HR_PIXCVec_008_511_219R_20240101T124346_20240101T124357_PIC0_01.nc
    xref_param_l2_hr_raster_file :
    SWOT_Param_L2_HR_Raster_20000101T000000_21000101T000000_20230817T100000_v302.rdf
    xref_reforbittrack_files :
    SWOT_RefOrbitTrackTileBoundary_Nom_20000101T000000_21000101T000000_20200617T193054_v101.txt, SWOT_RefOrbitTrack125mPass1_Nom_20000101T000000_21000101T000000_20200617T193054_v101.txt, SWOT_RefOrbitTrack125mPass2_Nom_20000101T000000_21000101T000000_20200617T193054_v101.txt
    utm_zone_num :
    11
    mgrs_latitude_band :
    S
    x_min :
    568300.0
    x_max :
    718600.0
    y_min :
    3900200.0
    y_max :
    4050500.0
    institution :
    CNES
    product_version :
    01
  • # Lisiting all available variables to select and visualize.
    variables = list(ds.data_vars)
    
    print("Variables within the NetCDF file:")
    for variable in variables:
        print(variable)
    Variables within the NetCDF file:
    crs
    longitude
    latitude
    wse
    wse_qual
    wse_qual_bitwise
    wse_uncert
    water_area
    water_area_qual
    water_area_qual_bitwise
    water_area_uncert
    water_frac
    water_frac_uncert
    sig0
    sig0_qual
    sig0_qual_bitwise
    sig0_uncert
    inc
    cross_track
    illumination_time
    illumination_time_tai
    n_wse_pix
    n_water_area_pix
    n_sig0_pix
    n_other_pix
    dark_frac
    ice_clim_flag
    ice_dyn_flag
    layover_impact
    sig0_cor_atmos_model
    height_cor_xover
    geoid
    solid_earth_tide
    load_tide_fes
    load_tide_got
    pole_tide
    model_dry_tropo_cor
    model_wet_tropo_cor
    iono_cor_gim_ka
    variable_name = 'wse'
    
    variable_to_plot = ds[variable_name]
    fig, ax = plt.subplots(figsize=(15, 10))
    contour = ax.contourf(variable_to_plot['x'], variable_to_plot['y'], variable_to_plot, cmap='viridis')
    cbar = plt.colorbar(contour)
    cbar.set_label('Water Surface Elevation (meters)', fontsize=12)
    ax.set_title('SWOT Raster 100m: Lake Mead Reservoir', fontsize=14)
    
    plt.show()

    raster_plot = ds.wse.hvplot.quadmesh(x='x', y='y', rasterize=True, title=f'SWOT Raster 100m: Lake Mead Reservoir')
    
    raster_plot.opts(width=1000, height=600, colorbar=True)

    Creating a Time Series

    SWOT Raster product does not include a time dimension, but it can be inserted by extracting from the file naming convention.

    1. Expand the time range of your earthaccess search to get an adequate range.
    2. Extracting the datetime from the s3 file naming convention then concatenating based on the new time dimension.
    3. Save the output as a new NetCDF file.
    time_results = earthaccess.search_data(
        short_name = 'SWOT_L2_HR_RASTER_2.0',
        bounding_box=(-114.502048,36.060175,-114.390983,36.210182),
        temporal =('2023-12-01 00:00:00', '2024-02-15 23:59:59'),
        granule_name = '*_100m_*',
        count =200
    )
    Granules found: 7
    fs_s3 = earthaccess.get_s3fs_session(results=time_results)
    raster = []
    for g in time_results:
        for l in earthaccess.results.DataGranule.data_links(g, access='direct'):
            raster.append(l)
    
    print(len(raster))
    raster
    7
    ['s3://podaac-swot-ops-cumulus-protected/SWOT_L2_HR_Raster_2.0/SWOT_L2_HR_Raster_100m_UTM11S_N_x_x_x_007_496_046F_20231211T024534_20231211T024555_PIC0_01.nc',
     's3://podaac-swot-ops-cumulus-protected/SWOT_L2_HR_Raster_2.0/SWOT_L2_HR_Raster_100m_UTM11S_N_x_x_x_008_205_109F_20231221T142037_20231221T142058_PIC0_01.nc',
     's3://podaac-swot-ops-cumulus-protected/SWOT_L2_HR_Raster_2.0/SWOT_L2_HR_Raster_100m_UTM11S_N_x_x_x_008_496_046F_20231231T233040_20231231T233101_PIC0_01.nc',
     's3://podaac-swot-ops-cumulus-protected/SWOT_L2_HR_Raster_2.0/SWOT_L2_HR_Raster_100m_UTM11S_N_x_x_x_009_205_109F_20240111T110542_20240111T110603_PIC0_01.nc',
     's3://podaac-swot-ops-cumulus-protected/SWOT_L2_HR_Raster_2.0/SWOT_L2_HR_Raster_100m_UTM11S_N_x_x_x_009_496_046F_20240121T201546_20240121T201607_PIC0_01.nc',
     's3://podaac-swot-ops-cumulus-protected/SWOT_L2_HR_Raster_2.0/SWOT_L2_HR_Raster_100m_UTM11S_N_x_x_x_010_205_109F_20240201T075048_20240201T075109_PIC0_01.nc',
     's3://podaac-swot-ops-cumulus-protected/SWOT_L2_HR_Raster_2.0/SWOT_L2_HR_Raster_100m_UTM11S_N_x_x_x_010_496_046F_20240211T170050_20240211T170111_PIC0_01.nc']
    def add_time_dimension(ds, file_path):
        # Extract filename from s3 file path
        file_name = file_path.split('/')[-1]
        # Extract date/time string from filename
        date_str = file_name.split('_')[-4][:15]
        # Convert the date string to a datetime object
        time_value = datetime.strptime(date_str, "%Y%m%dT%H%M%S")
        # Assign the time coordinate to the dataset
        ds.coords['time'] = time_value
        return ds
    datasets = []
    file_names = []
    for file_path in raster:
        with fs_s3.open(file_path, 'rb') as file:
            file_bytes = file.read()
        file_obj = io.BytesIO(file_bytes)
        dataset = xr.open_dataset(file_obj, engine='h5netcdf')
        dataset_with_time = add_time_dimension(dataset, file_path)
        datasets.append(dataset_with_time)
        file_names.append(file_path.split('/')[-1])
        dataset.close()
    # sorting the time dimension in correct order
    datasets.sort(key=lambda ds: ds.time.values)
    result_dataset = xr.concat(datasets, dim='time')
    result_dataset
    <xarray.Dataset>
    Dimensions:                  (x: 1549, y: 1549, time: 7)
    Coordinates:
      * x                        (x) float64 6.788e+05 6.789e+05 ... 8.336e+05
      * y                        (y) float64 3.9e+06 3.901e+06 ... 4.055e+06
      * time                     (time) datetime64[ns] 2023-12-11T02:45:34 ... 20...
    Data variables: (12/39)
        crs                      (time) object b'1' b'1' b'1' b'1' b'1' b'1' b'1'
        longitude                (time, y, x) float64 nan nan nan ... nan nan nan
        latitude                 (time, y, x) float64 nan nan nan ... nan nan nan
        wse                      (time, y, x) float32 nan nan nan ... nan nan nan
        wse_qual                 (time, y, x) float32 3.0 3.0 3.0 ... 3.0 3.0 3.0
        wse_qual_bitwise         (time, y, x) float64 8.053e+08 ... 8.053e+08
        ...                       ...
        load_tide_fes            (time, y, x) float32 nan nan nan ... nan nan nan
        load_tide_got            (time, y, x) float32 nan nan nan ... nan nan nan
        pole_tide                (time, y, x) float32 nan nan nan ... nan nan nan
        model_dry_tropo_cor      (time, y, x) float32 nan nan nan ... nan nan nan
        model_wet_tropo_cor      (time, y, x) float32 nan nan nan ... nan nan nan
        iono_cor_gim_ka          (time, y, x) float32 nan nan nan ... nan nan nan
    Attributes: (12/49)
        Conventions:                   CF-1.7
        title:                         Level 2 KaRIn High Rate Raster Data Product
        source:                        Ka-band radar interferometer
        history:                       2023-12-15T16:15:45Z : Creation
        platform:                      SWOT
        references:                    V1.1.1
        ...                            ...
        x_min:                         678800.0
        x_max:                         833600.0
        y_min:                         3900500.0
        y_max:                         4055300.0
        institution:                   CNES
        product_version:               01
    ds2 = result_dataset
    variable_name = 'wse'
    data = ds2[variable_name]
    
    fig, ax = plt.subplots(figsize=(10, 8))
    fig.set_tight_layout({'rect': [0.01, 0.01, 1.0, 1.0]})
    
    contour = ax.contourf(data.isel(time=0), cmap='viridis')
    cbar = plt.colorbar(contour)
    cbar.set_label('Water Surface Elevation (meters)', fontsize=14) 
    times = ds2.time.values
    
    # Function to update the plot for each time step
    def update(frame):
        ax.clear()
        contour = ax.contourf(data.isel(time=frame), cmap='viridis',)
        formatted_time = str(times[frame])[:-7]
        ax.set_title(f'Date: {formatted_time}')
        ax.set_xlabel('Longitude', fontsize=14)
        ax.set_ylabel('Latitude', fontsize=14)
        ax.text(0.5, 1.05, 'SWOT Raster 100M Lake Mead Reservoir', transform=ax.transAxes, ha='center', fontsize=14)
        return contour,
    
    # Creating a gif animation
    ani = animation.FuncAnimation(fig, update, repeat=True, frames=len(data['time']), blit=True, interval=3000)
    
    output = ('time_series.gif')
    ani.save(output, writer='pillow', fps=.5)
    
    with open(output,'rb') as f:
        display(Image(data=f.read(), format='gif'))
    
    plt.close(fig)
    ds2.close()
    <IPython.core.display.Image object>

    We can also save the time series to a new netcdf file.

    path = "SWOT_Raster_LakeMead_Time"
    os.mkdir(path)
    output_netcdf_path = os.path.join(path, "Output_Time.nc")
    result_dataset.to_netcdf(output_netcdf_path)
    
    print(f"Output complete: {output_netcdf_path}")
    Output complete: SWOT_Raster_LakeMead_Time/Output_Time.nc